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ABSTRACT 

We apply the joint lensing and dynamics code for the analysis of early-type galaxies, 
"cauldron", to a rotating N-body stellar system with dark matter halo which significantly 
violates the two major assumptions of the method, i.e. axial symmetry supported by a two- 
integral distribution function. The goal is to study how cauldron performs in an extreme 
case, and to determine which galaxy properties can still be robustly recovered. Three data 
sets, corresponding to orthogonal lines of sight, are generated from the N-body system and 
analysed with the identical procedure followed in the study of real lens galaxies, adopting an 
axisymmetric power-law total density distribution. We find that several global properties of 
the N-body system are recovered with remarkable accuracy, despite the fact that the adopted 
power-law model is too simple to account for the lack of symmetry of the true density distribu- 
tion. In particular, the logarithmic slope of the total density distribution is robustly recovered 
to within less than 10% (with the exception of the ill-constrained very inner regions), the 
inferred angle-averaged radial profile of the total mass closely follows the true distribution, 
and the dark matter fraction of the system (inside the effective radius) is correctly determined 
within ~ 10% of the total mass. Unless the line of sight direction is almost parallel to the total 
angular momentum vector of the system, reliably recovered quantities also include the angu- 
lar momentum, the V/cr ratio, and the anisotropy parameter 6. We conclude that the cauldron 
code can be safely and effectively applied to real early-type lens galaxies, providing reliable 
information also for systems that depart significantly from the method's assumptions. 

Key words: gravitational lensing — methods: N-body simulations — galaxies: elliptical and 
lenticular, cD — galaxies: kinematics and dynamics — galaxies: structure. 



1 INTRODUCTION 

Determining the structure of early-type galaxies and reliably pin- 
ning down their dark matter content is a crucial step in order to fully 
understand the formation and evolution processes of these systems. 

Within the currently favoured cosmological scenario, the 
ACDM paradigm, early-type galaxies are though t to be formed 
via hierarchical merging of lower mass galaxies {Toornrd 1 19771 : 
IWhite & Frerudll99ll : iBarnesIl 19921 : Icole et alj|2000h . While very 
successful in reproducing many observational features of ellipti - 
cal galaxies, including complex ones (see e.g. Ijesseit et al.ll2007I ). 
these formation models are still encountering difficulties in explain- 
ing the origin of the empirical scaling laws that correlate the globa l 
properties of early-type galaxies (see e.g. iRobertson et alj|2006h . 
Providing a reliable and detailed description of the mass density 
distribution, orbital structure and intrinsic properties of early-type 
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galaxies is therefore critical in order to enable stringent tests of 
galaxy formation models. 

For this reason, considerable effort has been devoted dur- 
ing the last decades towards the observation and the modeliza- 
tion of nearby early-type galaxies, by means of both stellar 
dynamics and X-ray studies, finding more or less strong evi- 
dence for a dark matter halo component (e . g. Fabbian c] 1 1989L 
Mould et al.l ll99C|, ISaglia. Berlin & Stiavellil 1 199 21 iBertin et all 



1994lFranx et alj| 19941 ICarollo et alJI 19911 Arnaboldi et al.lll996. 
RixetalJ 1 19971. iMatsushita et al.l 1 19981, [Loewenstein & White 



1999 LlGerhard et a l. 2001. Borriello et a l. 2003, R omanowskv et al.l 



20031 iHumphrev et al.l l200d iForbes et al.| |2008l and more ~re- 
cently the SAURON collaboration: see e.g. Ide Zeeuw et al 1 l2002l . 
lEmsellem et al.ll2004ICappellari et alj|2006t) . Both methods, how- 
ever, present some difficulties. In the case of stellar dynamics 
it is believed that some degeneracy can be present between the 
mass profile of the galaxy and the anisotropy of the stellar ve- 
locity dispersion tensor, which can be alleviated when higher or- 
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der velocity moments are available (see e.g. iGerhardl 19931) or by 
using physically motivated distribution functions (see e.g. iBertinl 
1200(1 for a discussion of this point). X-ray analyses, on the other 
hand, can seriously overestimate the total mass of the system if 
the assumption of hydrostatic equilibriu m for the hot gas does not 
hold, especially near the center (see e.g. IPellegrini & Ciottill2006l : 
ICiotti & Pellegrinil2008l) . 

A full understanding of the evolution of early-type galax- 
ies cannot be achieved without extending the study also to the 
mass density profile of objects at higher redshift (z > 0.1). This, 
however, has not been attempted until recently, due to obser- 
vational limitations and to the increased difficulty in extracting 
detailed kinematic information, which hinders traditional analy- 
ses based on stellar dynamics only. An effective solution in or- 
der to overcome these issues is constituted by a joint analy- 
sis which combines the constraints from stellar dynamics with 
the information obtained from gravitational lensing, when the 
early-type galaxy also happens to act as a lens with respect 
to a background source at hi g her redshift dTreu & Koopmansl 
20021: IKoopmans & TreJ 120021: iTreu & Koopmansl 120031 . |2004 



de Ven et alj|2008l) . IKoopmans et al .1 J2006h have successfully 



used this combined approach to analyse fifteen early-type lens 
galaxies (within a redshift range z = 0.06 - 0.33) disco vered in 
the Sloan Lens ACS Survey (SLACS, iBolton et alj[2006h as well 
as six systems between z ~ 0.5 and 1 from the Lenses Structure 
and Dynamics (LSD) Survey, showing that all of the examined sys- 
tems are well described by a power-law total density distribution 
very close to r~ 2 . The techniqu e for the joint lensing and dyn amics 
analysis has been expanded bv lBarnabe & Koopm ans (2007, here- 
after BK07) into a general and self-consistent method, completely 
embedded within the framework of Bayesian statistics, which puts 
constraints on the total density distribution of the lens galaxy by 
taking advantage of all the available data, i.e. not only the lensed 
image and a single stellar velocity dispersion measurement, but also 
the surface brightness distribution and the 2D kinematic maps (first 
and second projected velocity moments). 

Similar to other methods for the determination of the structure 
and internal dynamics of early-type galaxies, already mentioned 
above, the joint lensing and dynamics analysis also relies on a cer- 
tain number of assumptions. For ex ample, the simple and robust 
approach of IKoopmans et all J2006I) treats the gravitational lens- 
ing and the stellar dynamics as independent problems. The pro- 
jected mass distribution of the lens galaxy is modelled as a singular 
isothermal ellipsoid in order to determine the total mass within the 
Einstein radius, which is then used as a constraint for the dynami- 
cal model, where spherical symmetry and a specific prescription for 
the stellar orbital anisotropy are assumed. The more sophisticated 
framework of BK07 is designed to be very general and allows in 
principle for an arbitrary choice of the total potential. In practice, 
however, such freedom must balance against technical and compu- 
tational limitations. Therefore, in order to have a fast and efficient 
algorithm, the current implementation of the method, the cauldron 
codfl is restricted to axisymmetric potentials and two-integral stel- 
lar phase-space distribution functions (DFs). Under these hypothe- 
ses, it has been shown in BK07 that the method is capable of recov- 
ering with considerable accuracy the correct potential parameters 
and inclination angle, even in the presence of realistic noise. 

The point above raises the question of whether (and to what 
extent) the simplifying assumptions can be deemed valid for the as- 
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trophysical systems to which such methods are applied. In fact, real 
galaxies are not idealized objects, and there is no reason to expect 
them to be exactly axisymmetric (and neither triaxial ellipsoids) or 
to have two or three integrals of motion. Whereas axial symmetry 
generally seems to constitute a fairly good approximate descrip- 
tion for most early-type galaxies, a more detailed inspection, such 
as that allowed by the SAUR ON observations (see lEmsellem et all 
l2004lMcDermid et al .11 20061) , reveals a multitude of features indi- 
cating departure from axisymmetry, e.g. the presence of isophotal 
twist in the surface brightness distribution, minor axis rotation and 
kinematically decoupled cores. 

For this reason, in the present paper we apply our algorithm 
to the end-product of a two-component (stars plus dark matter) 
N-body simulation of a merger process, i.e. to a system which 
does not obey any restrictive prescription of symmetry, and there- 
fore violates the assumptions of the method. We aim to study how 
the cauldron algorithm performs when subjected to this kind of 
"crash-test", and which quantities can be robustly recovered even 
in such an extreme c ase. A similar approach has been followed by 
lThomasetai] l l2007l) , who have applied their three-integral axisym- 
metric Schwarzschild code to the study of non axisymmetric N- 
body merger remna nts, although without an y gravitational lensing 
information, and bv lMeneghetti et all t2007j) in the case of clusters 
of galaxies. 

The paper is organized as follows. In Section [2] we provide 
an overview of the cauldron algorithm for combined lensing and 
dynamics analysis. In Section [3] we summarize the properties of 
the N-body system that we use as lens galaxy. In Section|4]we de- 
scribe how the 2D maps of the simulated data with added realistic 
noise are obtained from the particle distribution. In Section [5] we 
apply the joint lensing and dynamics analysis to the simulated data 
and we present the results, which are then further discussed in Sec- 
tion[6] where also conclusions are drawn. 



2 THE CAULDRON ALGORITHM FOR JOINT LENSING 
AND DYNAMICS ANALYSIS 

In this Section we recall the main features of the cauldron algo- 
rithm. We refer the reader to BK07 for a fully detailed description 
of the method. 

The central tenet of a self-consistent joint analysis is to adopt 
a total gravitational potential (t> (or, equivalently, the total density 
profile p, from which O is calculated via the Poisson equation), 
parametrized by a set r\ of non-linear parameters, and use it simul- 
taneously for both the gravitational lensing and the stellar dynam- 
ics modelling of the data. As shown in BK07, these two modelling 
problems, while different from a physical point of view, can be ex- 
pressed in an analogous way as a single set of coupled (regularized) 
linear equations. For any given choice of the non-linear parame- 
ters, the equations can be solved (in a direct, non-iterative way) to 
simultaneously obtain as the best solution for the chosen potential 
model: (i) the unlensed source surface brightness distributions, and 
(ii) the weights of the elementary stellar dyn amics building blocks 
(e.g. orbits or two-integral c omponents, TICs, ISchwarzschildl 19791 : 
IVerolme & de Zeeuwl 20021) . This linear optimization scheme is 
consistently embedded in the framework of Bayesian statistics. As 
a consequence, it is possible to objectively assess the probability 
of each model by means of the evi dence merit function and, th ere- 
fore, to rank different models (see iMackavll 1991 1 19991 120031) . In 
this way, by maximizing the evidence, it is possible to recover the 
set of non-linear parameters r\ corresponding to the "best" potential 
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model. Here, in the context of Bayesian statistics, the "best model" 
means the most plausible model in an Occam's razor sense, given 
the data and the adopted form of the regularization (the optimal 
level of the regularization is also set by the evidence). 

Whereas the method is in principle extremely genera]^, its 
current practical implementation, which will be referred to as the 
cauldron algorithm, is more restricted in order to make it com- 
putationally efficient and applies specifically to axisymmetric po- 
tentials, 9(R, z), and two-integral DFs / = /(£, L-) (where E 
and L- are, respectively, the energy and the angular momentum 
along the rotation axis). Under these assumptions, the dynamical 
model can be constructed by making use of the fast BK07 numer- 
ical imple mentation of the two- integ ral Schwarzschild method de - 
veloped bv lCretton etakl d 19991) and lVerolme & de Zeeuwl J2002h . 
whose building blocks are not stellar orbits (as in the classical 
Schwarzschild method) but TICs The weights map of the optimal 
TIC superposition which best reproduces the observables is yielded 
as an outcome of the joint analysis. 

The cauldron algorithm has been s uccessfully te sted against 
the analytic power-law galaxy models of |E vans! J 1994i> . which re- 
spect by construction the assumptions of axisymmetry and two- 
integral DF, and afterwards has been employed fo r the detailed 
analy sis of the SLACS lens galaxy SDSS J2321-097 dCzoske et all 
I2008L hereafter C08). The latter is a case study which presents a 
benchmark data set particularly well suited to the needs of caul- 
dron, i.e. high-resolution HST/ACS images of the gravitational 
lensed source and of the surface brightness distribution of the lens 
galaxy, and 2D maps of the projected velocity moments of the lens 
galaxy derived from VLT-VIMOS observations. Therefore, the ob- 
servations of SDSS J232 1-097 will be used as a reference in or- 
der to generate the simulated observables for the present study (see 

Section|4j- 

It has been shown by the work of Koo pmans et alj d2006l) that 
a simple power-law model for the total density distribution pro- 
vides a satisfactory description for all of the SLACS lens galaxies 
examined so far. This has been further confirmed in the case of 
SDSS J232 1-097, where a fully self-consistent analysis was per- 
formed. In the present work, we aim to study the simulated galax- 
ies exactly as we would do for real objects, without assuming any 
a priori knowledge, and therefore we still adopt the same power- 
law model that has been used in C08. In particular, the total mass 
density distribution of the galaxy is taken to be a power-law strati- 
fied on axisymmetric homoeoids: 



p( m ) = fL < y < 3, 



(1) 



where po is a density scale, 7' will be referred to as the (logarithmic) 
slope of the density profile, and 
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where Co and ao are length-scales and q = co/ao. 

2 One could adopt for example a completely general pixelized potential for 
which the best profile is then determined via Bayesian statistics only by the 
data. 

3 A TIC can be visualized as an elementary toroidal system, completely 
specified by a particular choice of energy E and axial component of the an- 
gular momentum L-. TICs have simple i/R radial density distributions and 
analytic unprojected velocity moments, and by superposing them one can 
build f(E,L.) models for arbitrary spheroidal potentials (cf. ICretton et alj 
1 19991) : all these characteristics contribute to make TICs particularly valu- 
able and "inexpensive" building blocks when compared to orbits. 



The (inner) gravitational potential associat ed with a ho- 
moeoi dal density distribution p(m) is given by IChandrasekhaJ 
l l 19691) formula. In our case, for 7' i= 2, one has 
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There are three free non-linear parameters in the potential to 
be determined via the evidence maximization: O ( or equivalently, 
through equation [B4] of BK07, the lens strength a ), the slope 7' 
and the axial ratio q. When required by the data, it is straightfor- 
ward to include a core radius R s in the density distribution. Beyond 
the previously mentioned parameters, there are four additional pa- 
rameters which determine the geometry of the observed system: 
the position angle j? pa , the inclination i and the coordinates of the 
centre of the lens galaxy with respect to the sky grid. The position 
angle and the lens center can usually be accurately determined by 
means of a preliminary exploration and kept fixed afterwards in or- 
der to reduce the number of free parameters. Finally, for a proper 
modelling of the lensed image it can be necessary to include two 
more parameters (shear strength f and shear angle j? f ) in order to 
account for external shear. 

A curvature regularization (as described in ISuvuetal.ll2006l 
and appendix A of BK07) is adopted for both the gravitational 
lensing and the stellar dynamics reconstructions. As discussed in 
BK07, the initial guess values of the hyperparameters (defining the 
level of the regularization) are chosen to be quite large, since the 
convergence to the maximum is faster when starting from an over- 
regularized system. 



3 THE N-BODY SYSTEM 

The model galaxy that we use as lens in the present work is 
the end-product of a numerical N-body simulation of a dissi- 
pationless merging between two equal-mass spherical galaxies 
embedded in their dark matter halos. The simulation, run with 
the treecode fvfps (Fortran Version of a Fast Poisson Solver ; 
lLondrillo. Nipoti & Ciottill2003trNipoti. Londrillo & Ciotti|[2003h . 
has been already presented and described in a previous paper 
l lNipoti. Londrillo & Ciottil2007[) . in which it is named E25o. Here 
we recall the main propertie s of the end-product of this simulation, 
while we refer the reader to iNipoti et al j d2007l) for a detailed de- 
scription of the initial conditions. 

The simulation end-product is a nearly ellipsoidal virialized 
system, comprising a stellar component and a dark-matter compo- 
nent. The total number of particles is ~ 1.2 x 10 6 , and all particles 
have the same mass. The stellar component has total mass ~ 2 M t 
and angle-averaged half-mass radius ~ 3.8 r», where M, and r, 
are the code length and mass unit. The dark matter component is 
more massive and more extended than the stellar component, hav- 
ing mass ~ 10 M, and half mass radius ~ 19.3 r„. The system has 
a virial velocity dispersion ~ 0.55(GM,/r t ) l/2 and non-vanishing 
total angular momentum L, corresponding to a value A ~ 0.07 in 
terms of the spin parameter A = |£ t otl 1 ' 2 l|L||G" 1 M^, , where E lot is 
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Figure 1. Axis ratios b/a (solid curves) and c/a (dotted curves) as func- 
tions of radius for the stellar (thin curves) and total (thick curves) density 
distributions of the N-body system used as lens (we assumed r„ = 1.969 
arcsec). 

the total energy, and M tot ~ 12 M, is the total mass. Of course, the 
model can be rescaled to represent physical systems of any size and 
mass by choosing proper values of M, and r„. 

We define the minor- to-major (c/a) and intermediate-to-major 
(b/a) axis ratios of the system at a radius r as the corresponding 
axis ratios of the inertia ellipsoid of part icles within an ellipsoid of 
angle - averaged radius r (for details, see lNipoti. Londrillo & CiottH 
l2002h . We find that the system is triaxial, with direction of the prin- 
cipal axes and axis ratios depending on radius. Figure [T] shows b/a 
(solid curve) and c/a (dotted curve) as functions of radius for the 
stellar (red) and total (black) density distributions within about the 
half-mass radius of the stellar component (we assumed r, = 1.969 
arcsec, see Section[4]l. Both the stellar and the total distribution are 
strongly triaxial at r > 5 arcsec and mildly triaxial (almost prolate) 
at r < 5 arcsec. The angle- averaged total density and mass pro- 
files of the model are plotted as solid black curves in Figs.[3]and|4] 
respectively. 



4 CONSTRUCTION OF THE OBSERVABLES 

In this Section we detail how the simulated observables (or "mock 
data") are generated from the dark matter and stellar particle distri- 
bution taken from the virialized end-product of the N-body simula- 
tion described in Section[3] 

In order to convert the N-body simulation in a realistic data set 
with plausible physical characteristics (effective radius, redshift of 
lens galaxy and source, Einstein radius), we use as a reference the 
actual lens galaxy SDSS J2321-097, for which both data analysis 
and joint lensing and dynamics study are presented in C08. There- 
fore, we will adopt for the noise level and the sampling of the data 
(i.e. size and binning of the data grids) the corresponding values of 
SDSS 12321-097. 

First, we impose for the simulated galaxy the same redshift of 
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Figure 2. Mock source to be lensed. This is an elliptical Gaussian brightness 
distribution with <r x = 0.2 arcsec, cr lJ = 0.1 arcsec, position angle #pa = 
115° and with a slight offset with respect to the lens center. 

SDSS 12321-097, i.e. a = 0.0819. The adopted source (Fig. [2j to 
produce the artificial lensed image is an elliptical Gaussian distri- 
bution, slightly offset with respect to the center of mass of the lens 
galaxy and (again in analogy with SDSS 1232 1-097) located at a 
redshift z s = 0.5342. We then fix the values r, = 3 kpc (^ 1.969 
arcsec at redshift zi) and M* = 3 X 10 M e for, respectively, the 
length and mass units of the simulated galaxy. This setup produces 
a galaxy that, when is observed along an arbitrary line of sight, dis- 
plays realistic values for both the effective radius R,, (~ 5-6 arcsec, 
corresponding to ~ 8 kpc at the redshift z\), and for the Einstein ra- 
dius ^Ei nst (~ 2 arcsec). 

With this choice of r, and M,, the simulated galaxy is a mas- 
sive system of total mass 3.6 x 10 12 Mq, with a dark halo which 
extends up to several hundred kpc. In analogy with the data set of 
SDSS 1232 1-097, the spatial coverage of our data is confined to 
the very inner regions of the galaxy, approximately corresponding 
to /? c /2. However, a fair amount of information comes also from 
more distant regions of the system which are seen in projection 
along the line-of-sight (see discussion in C08). 

We define for the galaxy an orthogonal reference frame with 
the origin in the center of mass of the simulated system; the z axis 
is oriented along the direction of the total angular momentum of the 
stellar component. Since all the observables are quantities projected 
in the sky plane, we select three orthogonal lines of sight along 
which the simulated galaxy is assumed to be observed. In this way, 
from a single simulation, we obtain three different data sets. The 
first line of sight is chosen to be approximately oriented along the 
z axis of the galaxy, and therefore in the following we will refer 
to the corresponding projection as the .w/-plane projection or, for 
simplicity, the "face-on" projection. The remaining two projections 
are taken along mutually orthogonal axes both perpendicular to the 
first line of sight and will be called the yz- and zx-plane projections 
or, for brevity sake, the "edge-on" projections. 

The circularized half-light radii for the three projections are 
found to be R e = 5.227 arcsec (i/z-plane), 5.359 arcsec (zjc-plane) 
and 6.156 arcsec (xj/-plane). 

When considering any one of the three projections, we will 
always indicate as z' the direction of the particular line-of-sight. 
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For each projection, the observables are then calculated following 
the same procedure: 

(i) Lensed image. All the particles (both stellar and dark mat- 
ter ones), properly weighted by the respective masses, are cast in 
projection along the line-of-sight and binned in a 2D-grid along 
the sky plane in order to generate from the total density p m a pro- 
jected density map E 1o1 . The latter, once normalized to the critical 
density Z cl = (c 2 /4nG) CD s /AjAis), a quantity which depends only 
on geometry of the systetrQ, yields the convergence field k on the 
projection plane. The two-dimensional Poisson equation V 2 i/e = 2k 
relates the convergence to the projected potential if/, whose gradi- 
ent im mediately gives the deflection ang le vector field: a = Wifi 
(see e.g. lSchneider. Ehlers & Falcolll992 ). 

We make use of a square grid of 3000 x 3000 bins for the con- 
vergence. The grid size (equivalent to ~ 40 arcsec) is such that 
it contains, in projection, about half of the total number of parti- 
cles. Including more distant particles does not have any discernible 
effect on the outcome, except slowing down all the related calcula- 
tions, since the number of bins must be increased in order to keep 
the resolution constant. The Poisson equation is then solved via 
fast Fourier transform, using the freely available package FFTW 
dFrigo & Joh nson 2005), and enabling us to obtain, from the con- 
vergence (appropriately padded in order to avoid numerical issues 
with the Fourier transform), the two components of deflection an- 
gle a. With these deflection angle maps and the mock source de- 
scribed above, the cauldron code can straightforwardly generate 
the lensed image, already convolved with the HST/ ACS F814W 
point spread function (PSF) obtained with tiny tem dKristll 19931) . 
The lensed image is constructed (in analogy with the lensing data 
for SDSSJ2321-097) on 100 x 100 grid, with each pixel corre- 
sponding to 0.05 arcsec. Finally, a noise distribution based on the 
covariance maps of the HST images of SDSS J2321-097 is added 
to the lensed image map in order to produce the final data set. 

The obtained lensing data sets for the three projections are shown 
in the upper right panel of Figs[5][7]and[9] 

(ii) Surface brightness distribution. We assume that the stellar 
mass-to-light ratio is independent of position. The stellar particles 
are cast along the chosen line-of-sight on a 2D-grid in the sky plane. 
Such grid is padded and oversampled by a factor of 3 with re- 
spect to the final grid adopted for this observable (see later). This 
is necessary since, in order to take into account the effect of seeing, 
we have to convolve this quantity with the PSF (here modelized 
as a Gaussian distribution of FWHM = 0.10 arcsec). The map is 
then resampled, generating the surface brightness distribution of 
the galaxy on a 50 x 50 grid (1 pixel = 0.10 arcsec). Mock noise 
is then added consistently with the corresponding variance maps 
of SDSS J232 1-097. These images are to first-order equivalent to 
HST-NICMOS images. 

The obtained data sets for the surface brightness distribution of 
the three projections are shown in the upper left panel of Figs [6] [8] 
and[H 

(iii) Line-of-sight projected velocity moments. The procedure to 
generate the kinematic maps is similar to the one described above 
for the surface brightness distribution, including the oversampling 
and convolution with the Gaussian PSF (but with a broader FWHM 
= 0.90 arcsec, typical for ground-based observations with the 
VLT). 



4 D s , D\ and D[ s are the angular diameter distances from the observer to 
the source, from the observer to the lens and from the lens to the source, 
respectively. 



Table 1. Recovered non-linear parameters of the best power-law models for 
the three data sets (obtained as projection of the simulated systems on the 
three orthogonal planes indicated in the columns): inclination r' (in degrees), 
shear strength £ and angle ff( (in degrees), lens strength ag, logarithmic 
slope Y and flattening q. 
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The only difference is that now the velocities projected along the 
line-of-sight (v z > and u 2 , for the two kinematic maps respectively) 
associated with each particle are summed up on each cell, produc- 
ing the unweighted maps to be used by cauldron. These maps can 
be divided by the surface brightness distribution sampled on the 
same grid (the "kinematic" grid) in order to obtain the weighted 
maps, i.e. the quantities (ty > and (u 2 ). The projected velocity dis- 
persion is obtained as <x 2 = {v 2 _,) - {v-j) 2 . 

Due to the challenges of spectroscopic observations of distant 
early-type galaxies, the kinematic grids of SDSS J232 1-097 have 
only 9x9 elements, with 1 pixel = 0.67 arcsec (i.e. the VLT- 
VIMOS IFU fiber size). Since we are mimicking those observation, 
we adopt the same grid. As usual, the mock noise is added accord- 
ing to what we know from SDSS J2321-097. 

The obtained kinematic data sets for the three projections are 
presented in the upper central (for the line-of-sight velocity) and 
upper right (for the line-of-sight velocity dispersion) panels of 
FigsHHUandQ! 



5 ANALYSIS AND RESULTS 

In this section we illustrate the results of the joint lensing and dy- 
namics analysis performed with cauldron on the three data sets 
generated (as described in Section [4) from orthogonal projections 
of the simulated galaxy. 



5.1 Recovered structure 

As discussed in Section [2] we have adopted an axisymmetric 
power-law profile as a model for the total density distribution. The 
recovered non-linear parameters of the best models for the three 
data sets are reported in Table [T] The non-linear parameters listed 
in the Table are the inclination ( (in degrees), the shear strength £ 
and angle j?^ (in degrees), the lens strength a , the logarithmic slope 
Y and the flattening q (a q larger than 1 denotes a prolate axisym- 
metric shape). In order to speed up the optimization routine, the 
best values for the lens center and galaxy position angle were de- 
termined via preliminary runs, and afterwards were kept fixed (this 
procedure is commonly used in this kind of analysis: see e.g. C08). 
The core radius R s was initially set free as an additional non-linear 
parameter, but, in order to make the optimization faster, it was later 
fixed to a negligibly small value after verifying that the introduc- 
tion of this new parameter did not lead to an increase of the value of 
the evidence merit function (and therefore in a Bayesian sense the 
additional complexity was not justified). The optimization routine 
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Figure 3. Angle-averaged total density profile of the simulated system 
(black line) compared with the density profiles of the recovered best mod- 
els (see Table[TJ for the three orthogonal projections data sets. The i/z-plane 
model (red line) has a logarithmic density slope Y = 2.214, the z;t-plane 
model (blue line) has a slope yt x = 2.078, and the Ary-plane model (green 
line) has a slope y' nJ = 2.234. 
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Figure 4. Angle-averaged total mass distribution of the simulated system 
(black line) compared with the total mass profiles of the recovered best 
models (see Table[T} for the three orthogonal projections data sets: f/z-plane 
model (red line), z-t-plane model (blue line) and .vy-plane model (green 
line). 



yields also the best value for the three hyperparameters which set 
the ideal level of regularization. 

The reconstructed observables for gravitational lensing and 
stellar dynamics that correspond to the best model for the yz-plane 
data set are presented and compared to the data in Figs [5] and [6] 
respectively. Analogously, Figs. |7|8| and |9|10l show the same quan- 
tities for the zx-plane and xj/-plane (i.e. face-on projection) data 
sets. 

It is apparent that the residuals of the reconstructed lensed 
image are fairly large, in particular in the cases of the zx- and 
xi/-plane projections, where moreover the reconstructed source ap- 
pears patchy and unrealistically pixelized, despite the regulariza- 
tion. The same effects, while less pronounced, are discernible also 
in the case of the ^z-plane projection, where the reconstruction was 
most successful. This is a clear indication that the underlying total 
density distribution is actually more complex or inhomogeneous 
than the simple power-law profile that we are using as a model (see 
Section 1531 for a more extended discussion) and that the model is 
unable to de-lens all lensed images simultaneously into a single 
well-defined source. 

The surface brightness distribution and the kinematics appear 
to be reasonably well reconstructed. However, in the inner region 
the reconstructed velocity dispersion is more peaked than in the 
data, with the possible exception of the zx-plane data set (the one 
for which the recovered slope is the most shallow, cf. Table [TJ. To- 
gether with the faint central image visible in all the lensing data 
sets, this indicates that the total density distribution of the system 
is shallower in the central regions than in the outer regions, as con- 
firmed by the direct analysis of the angle-averaged density profile 
of the object (see Fig. [3}. 

Despite several indications from the reconstructed observables 
that the adopted power-law model is probably oversimplified for 



the data at hand, we find that such a model still provides a satisfac- 
tory description of the essential features of the system, and various 
important physical quantities are robustly recovered. 

The angle-averaged total density profiles corresponding to the 
best models of Table [T]for the three data sets are plotted in Fig. [3] 
and compared with the true profile of the simulated system. The 
density slope, which is very close to y' ~ 2.2, is quite accu- 
rately recovered beyond the inner ~ 0.5 arcsec. The zx-plane model 
presents a slightly shallower profile (y' ~ 2.1), although the dis- 
crepancy is small. The ;q/-plane model, which in general provides 
the worst recovery of the true quantities (as it will be further dis- 
cussed in this Section), has a density normalization which is lower 
than the real one, but manages to catch almost perfectly the cor- 
rect density slope. In the inner half arcsecond, the density profile 
of the simulated galaxies becomes shallower (p ~ 1/r), a feature 
that the power-law model is obviously unable to capture, although, 
as already discussed, we have other signals (e.g. from the velocity 
dispersion) that the model breaks down in that approximate region, 
and that this is due to the presence of a break in the density pro- 
file. Slightly more sophisticated models, such as a double power- 
law or a single power-law with an added break radius, have also 
been explored, but the reconstruction of the observables does not 
improve significantly and the additional complexity is therefore pe- 
nalized by the Bayesian evidence in favour of the single power-law. 
In light of what we know from the direct examination of the simu- 
lated object, the reason for this failure is that the transition from the 
p ~ 1/r 2 region to the p ~ 1/r region is too abrupt to be adequately 
reproduced by these models, and therefore they do not perform sig- 
nificantly better than the simpler single power-law. 

Whereas the density slope is a substantially unharmed sur- 
vivor of the crash test, the recovered flattening and inclination angle 
are not reliable parameters in case the systems deviate too dras- 
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Table 2. Dark matter fraction within a sphere of radius r = R c (first column) 
and within the line-of-sight oriented cylinder of radius R = R e (second 
column) for the three best models. The corresponding quantities for the true 
system within the same radius are also presented. 





DM fraction 


DM fraction 




(sphere) 


(cylinder) 


i/z-plane model (R c = 5'.' 23) 


0.20 


0.29 


true system (at the same radius) 


0.16 


0.33 


z.v-plane model (R e = 5'.'36) 


0.24 


0.37 


true system (at the same radius) 


0.16 


0.34 


.ty-plane model (if c = 6'.' 16) 


0.19 


0.25 


true system (at the same radius) 


0.20 


0.35 



tically from the model's assumptions. One should note, however, 
that these quantities are only properly defined for an axisymmetric 
object, and therefore do not have a straightforward interpretation 
when directly applied to a simulated system which is approximately 
triaxial and whose axis ratios also change as a function of radius (as 
seen in Section[3j- The best models for the i/z-plane and xy-plane 
data set both give an inclination close to ~ 55°, so there is no sign 
that the latter is interpreted as a face-on system. The zx-plane best 
model, on the other hand, turns out quite correctly to be an almost 
edge-on system (;' =s 90°). As for the flattening, it appears that the 
axisymmetric model, faced with the insurmountable problem of tri- 
axiality, tends to adopt a density profile close to spherical (almost 
exactly spherical in the case of the z/z-plane data set, slightly prolate 
or oblate in the cases of the zx- and xi/-plane data sets, respectively). 



5.2 Total mass distribution and dark matter fraction 

Closely connected to the density profile is the (angle-averaged) to- 
tal mass distribution, plotted in Fig.[4]for the three best models and 
the true system. We find that in the case of the edge-on projections 
models (yz- and zx-planes) the mass profile is very well reproduced 
within a few percent. The mass profile of the face-on projection 
model is instead underestimated of about 25 percent within a sphere 
of radius ~ 5 arcsec: this is a consequence of the too low density 
normalization which is recovered for this model, as previously seen 
in Fig. [3] Because of the tight constraints imposed by the lensing, 
the total projected mass enclosed within /?Ej nst is within a few per- 
cent from the correct value for all three models. 

A quantity of extreme interest in the study of early-type galax- 
ies is the dark matter fraction of these objects. It is therefore impor- 
tant to be able to test how reliable is the cauldron method in esti- 
mating this parameter also for systems that defy its assumptions of 
axisymmetry and two-integral DF. Since the method only provides 
a total density distribution, however, it is necessary to make fur- 
ther assumptions to be able to constrain the stellar density profile. 
In order to limit as much as possible the arbitrariness of such as- 
sumptions we adopt, as in the analysis of real galaxies, the so-called 
"maximum bulge" approach (cf. C08). This consists in maximizing 
the contribution of the luminous component (which is obtained as 
an output of the best model reconstruction), i.e. maximally rescal- 
ing the stellar density distribution without exceeding the total den- 
sity distribution, assuming that the stellar mass-to-light ratio is in- 
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Figure 5. Best model lens image reconstruction for the i/z-plane data set 
(the generation of the simulated observables is detailed in SectionfJ}. From 
the top left-hand to bottom right-hand panel: reconstructed source model; 
simulated noisy data showing the lensed image; lensed image reconstruc- 
tion; residuals. 
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Figure 6. Best dynamical model for the i/z-plane data set. First row: sim- 
ulated noisy surface brightness distribution, projected line-of-sight veloc- 
ity and line-of-sight velocity dispersion. Second row: corresponding recon- 
structed quantities for the best model. Third row: residuals. 



dependent of position. This method gives a lower limit for the dark 
matter fraction, provided that the model's assumptions hold true. 

Under the maximum bulge hypothesis, we study both the vol- 
ume mass ratio (i.e. the dark matter fraction within a sphere of ra- 
dius taken to be equal to the effective radius of the considered data 
set) and the projected mass ratio (i.e. the dark matter fraction within 
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Figure 7. Same as Fig. [5] but relative to the jjr-plane data set. 



Figure 9. Same as Fig. [5] but relative to the ja/-plane data set. 
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Figure 8. Same as Fig. [6] but relative to the zjc-plane data set. 



Figure 10. Same as Fig. [6] but relative to the jcy-plane data set. 



a cylinder oriented along the line-of-sight with a radius equal to 
R c ). The results of this analysis for the three best models are sum- 
marized in Table [2] which also shows the corresponding quantities 
for the true system (at the appropriate radii). We find that the dark 
matter fraction is remarkably well recovered for all three models, 
and it is within 10% (in total mass) of the correct value for both the 
volume and the projected mass ratio. The largest discrepancy (10% 
within the cylinder) is found for the face-on model, as usual the 
most problematic one, while the j/z-plane model (which is the one 
that best reproduces the observables, in particular the lensed im- 
age) manages to accurately recover the dark matter fraction within 
< 5% of the correct value. 



It must be noted that we assumed a position-independent stel- 
lar mass-to-light ratio both in constructing the surface-brightness 
map from the N-body system and in estimating the maximum stel- 
lar mass for the best-fit models. In real galaxies the stellar mass- 
to-light ratio might depend on position, though the effect of a non- 
uniform stellar mass-to-light ratio is not expected to be stro ng based 
on observed colour gradients (e.g. lKronawitter et al.ll200ol) . 

5.3 Reconstructed source 

As it can be immediately seen by an examination of Figs|5]|7]and|9] 
the lensed image reconstruction for the three models is far from the 
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noise level, and the reconstructed source - with the partial excep- 
tion of the j/z-plane model - has little in common with the simulated 
source (Fig.O used to generate the lensed images. Even assuming 
that one has no information about the actual source (as it would 
be the case for real data sets) it is evident that the reconstructed 
sources are unrealistically irregular and pixelized. We remark that 
this is a consequence of having an underlying density distribution 
for the system which is more complex and less homogeneous than 
the assumed power-law density model. As the present test neatly 
displays, gravitational lensing is very sensitive to the features of 
the potential within the Einstein radius, and even small inhomo- 
geneities and departures from the assumed model will generally 
have a detectable effecjf]. Significantly, such large residuals have 
never been encountered so far when analysing real lens galaxies, 
with all the examined SLACS lenses being accurately reconstructed 
by means of single p ower-law models, poss i bly with the inclu- 
sion of external shear jKoopmans et aljfeood , ICzoske et alj|2008l 
Barnabe et al., in preparation). This strengthens the statement that 
the simulated system under analysis, while not being utterly un- 
realistic, constitutes an extreme case and is therefore well suited 
for the "crash test" of the code that we aim to conduct. It should 
be emphasized that, if these were real data sets, from the mere vi- 
sual inspection of the results of the best lensed image reconstruc- 
tion provided by cauldron, we would already be able to conclude 
that we are dealing with a unusually complex galaxy, which would 
surely deserve a more sophisticated modelling once the preliminary 
study is completed. 

In order to better understand the limitations of the axisym- 
metric approach and the cause for the large residuals in the recon- 
structed image, the same data set (i.e. i/z-plane projection) has been 
analysed by means of the adaptive Bayesian strong lensing code 
of VK08, which can account for small corrections in the projected 
potential, and therefore for departures from symmetry and for the 
presence of dumpiness and substructure in the convergence, but 
which does not include any constraint from the dynamics. Starting 
from the best model axisymmetric potential of Table Q] and then 
introducing small potential corrections and letting the parameters 
vary, it is found that the lensed image can be reconstructed down 
to the noise level, and the source is quite accurately recovered, by 
making use of an adaptive grid, as shown in Fig. QT] As it is vis- 
ible in the convergence map of Fig. QTJ (bottom-right panel), there 
are two main overdensities (with respect to the best model density 
profile), located on opposite sides with respect to the center. These 
overdensities are at least an order of magnitude too weak to be in- 
dicative of the presence of a genuine localized massive dishomo- 
geneity in the simulated system (cf. VK08 and in particular their 
Figs. 8 and 9), which is rather smooth and does not present massive 
substructures. These features in the convergence map are therefore 
due, rather than to dumpiness, to the deviation of the true (pro- 
jected) mass distribution from the simple assumption of elliptical 
shape. 

5.4 Recovered dynamical quantities 

A reliable knowledge of the dynamical structure of early-type 
galaxies would constitute a valuable asset for all the formation and 

5 As a consequence, gravitational lensing can actually be used to quan- 
tify the level of mass substructure in massive galaxie s, through their effect 
on hi ghly-magnified arcs and Einstein rings (see e.g. IVegetti & Koopmanj 
I2008L hereafter VK08). 



Table 3. Recovered dynamical quantities for the three best models (last 
three columns) compared with the true values directly calculated from the 
N-body system (second column). The dynamical quantities are L. (in units 
of kpc km s - '), the V/cr ratio, and the three global anisotropy parameters 6, 
fi and y. See text for a more exhaustive description. 





true value 


i/z-plane 


zx-plane 


jcy-plane 




-143.1 


-142.6 


-143.3 


16.8 


V/cr 


0.170 


0.152 


0.214 


0.130 




0.301 


= 


= 


= 


7 


0.208 


-0.470 


-0.645 


-0.986 


6 


0.219 


0.190 


0.244 


0.330 



evolution models. Therefore, in this Section we investigate how ac- 
curately the axisymmetric cauldron code can recover the essential 
dynamical characteristics of the simulated system. We expect com- 
parable performances (and typically better ones) when the code is 
applied to real galaxies. 

Since the dynamical modelling block of cauldron (see Sec- 
tion|2j does not employ an actual orbit-superposition method, a de- 
tailed analysis of the orbital families of the galaxy is beyond reach. 
However, we are able to study fundamental global dynamical quan- 
tities such as angular momentum, V/cr and anisotropy parameters. 
The results of this study for the three models, compared with the 
true quantities calculated directly from the simulated system, are 
listed in TablefJ] 

The total angular momentum along the principal rotation axis, 
L-, has been calculated for the three models and the simulated sys- 
tem within the same cylindrical region of radius R = 5" and height 
z = 5" (this is a square box in the meridional plane, of linear size 
comparable to the effective radius; all the dynamical quantities have 
been calculated within this region). As seen in Table[3] L- is accu- 
rately recovered for the two edge-on projection models, with dis- 
crepancies of less than 10% from the correct value, despite the fact 
that the line-of-sight velocity maps are considerably nois)[j§ On the 
contrary, in the case of the face-on projection model, where little 
or no information about the rotation is available from the kinematic 
data (see top-middle panel of Fig. [10] the map displays essentially 
no rotation) the recovered L z is obviously incorrect, being of the op- 
posite sign and, more importantly, quite close to zero. This clearly 
highlights the great importance of the information enclosed in the 
kinematic maps, despite the fact that such maps are usually much 
more noisy and considerably more coarsely sampled, in terms of 
number of pixels, than the surface brightness or lensed image maps. 
On the opposite side, this result also cautions us in being aware of 
the possible shortfalls of the method when studying galaxies whose 

data show no discernible rotation. 

The global quantity V/cr (see lBiniievll2005l : ICappellari et al .1 

120071) is an indicator of the importance of rotation with respect 
to the random motion. We find here that a value not too far from 
the correct one is recovered for the edge-on projection models, al- 
though the discrepancies are larger than in the case of the angular 
momentum (the difference is of order 10% for the i/z-plane model 



The mock galaxy displays also some rotation along the orthogonal axes, 
although the angular momentum along these directions is only of order 
< 1% of L-. Because of its intrinsic properties, the model is not capable 
of reproducing this kind of rotation (in other words, L x and L y are always 
by construction). 
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Figure 11. Result of the linear source and potential reconstruction for the i/z-plane dataset, as obtained by applying the adaptive Bayesian lensing code of 
VK08. The first row shows, from left to right, the lensed image data set, the reconstructed image and best source. The second row presents, from left to right, 
the image residuals, the total potential correction and the substructure convergence. 



and 25% for the zx-plane model), and even the face-on projection 
model, while underestimating the importance of rotation of about 
25%, gives a quite reasonable result. 

The anisotropy distribution is another very relevant dynami- 
cal quantity, often considered an important indicator of the galaxy 
formation processes. Since the simulated system is severely non- 
spherical, however, the anisotropy profile cannot be reliably de- 
scribed and compared to the models by making use of a sim- 
ple radial parameter such as the commonly used f} r (r) = 1 - 
°"?an('")/' T i 2 ad( r )' wnere °"ian an d °"iad are the tangential and radial ve- 
locity dispersioni, respectively. Instead, a more robust and consis- 
tent indic ator is provided by the th ree global anisotropy parameter s 
defined in lCappellari et al] d2007h and lBinnev & Tremaind d2008l) : 



n- 

= 1 _ ^ w 



2n- 



n«« + tl 



2j3-r 

2-7 ' 



where 



(6) 
(7) 
(8) 

(9) 



and cr k denotes the velocity dispersion along the direction k at any 
given location in the galaxy. For a two-integral DF a\ = a 1 , every- 



where, which implies Tl RR = IT.-, so that the value of p is always 
zero. Therefore, due to the simplifying assumptions on the distri- 
bution function, our method will generally fail in recovering the 
true p value when analysing a more complex system. We clearly 
observe this in the present study, where (3 = 0.301 for the simu- 
lated system (see again Table [3}- Not too unexpectedly, since the 
global anisotropy parameters are related, this has disrupting con- 
sequences also on the recovered y, which is always negative, indi- 
cating mild tangential anisotropy in the models, while the mildly 
radially anisotropic N-body system has a positive y. This discrep- 
ancy is confirmed by inspecting the radial behaviour of a kinetic 
energy proxy for/3 r O), i.e. the quantity fi K (r) = 1 - K tm (r)/K T!li (r), 
where A" tan and K mlj are, respectively, the spherically-averaged tan- 
gential and radial components of the kinetic energy. Interestingly, 
however, the global parameter 5 is remarkably robust, particularly 
for the two edge-on projection models, with a discrepancy of order 
~ 15% from the true value. 



5.5 Uncertainties 

In this section we present the statistical uncertainties on the model 
parameters calculated by making use of the identical procedure 
which is used when dealing with real systems. 

Bayesian statistics represents a powerful tool for data analysis, 
model comparison and model parameters constraining. A major im- 
provement in this direction has been made with the introduction of 
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Figure 12. Posterior probability distribution for the non-linear parameters of the best power-law model for the yz-plane projection data set, as obtained from 
the Nested Sampling evidence exploration. The width of the posterior probability distribution has been calculated for each of the parameters, by considering 
the region around the peak which contains 99% of the probability, yielding: i = [57.81,58.45], £ = [0.05273,0.05455], & ( = [-26.55,-25.47], a = 
[0.6512,0.6550], / = [2.231,2.244], q = [1.008, 1.023]. 



the Nested Sampling techn ique d eveloped by Skillind 1 2004 ) (see 



also Sivia & Skillind2006j and e.g. lMukheriee. Parkinson & Liddld 
12006 for an astrophysical application). This method provides a 
computationally efficient way to calculate the total marginalised 
evidence, which is the key quantity for model comparison, and in 
addition yields other valuable information, such as the posterior 
probability density distributions, the mean values and the standard 
deviations for each of the model parameters. The Nested Sampling 
technique has been applied in the context of lensing and model- 
comparison by VK08. 

Bayesian statistics requires to formalize o ne's assumption s by 
defining priors P on the parameters 77; (see e.g. iMackavll 1 9921) . We 
choose the priors to be uniform in a symmetric interval of size Sr/, 
around the best recovered values r] bi , that is: 



constant for |;/ b i - 77; | < 6t]i 



(10) 







for l'7b,i - ^il > 



In order to make sure that the priors include the bulk of the evi- 
dence likelihood, very conservative estimates of the intervals Sr]i 
are obtained by means of fast preliminary runs. 

The posterior probability distributions (PPDs) for each pa- 
rameter obtained from the Nested Sampling analysis are shown in 
Fig-E]for one of the models (i/z-plane projection). Within the con- 
text of Bayesian statistics, each PPD histogram quantifies the error 
for the considered parameter given the data and all the assumptions 



(i.e. under the hypothesis that the adopted model is the correct one); 
since the many pixels in the data sets provide a lot of constraints, 
these errors are typically very small, as is the case in the plot pre- 
sented here. It should be noted that, because of the marginalization 
over all the parameters except the one under analysis, the PPD usu- 
ally provides the most conservative estimate of statistical errors. At 
the same time, however, due to projection effect arising from the 
marginalization itself, the mean value of the parameter obtained 
from the PPD can be significantly skewed with respect to the best 
recovered value of the corresponding parameter obtained from the 
best model (P. Marshall, private communication). 

However, the statistical errors cannot take into account or give 
an estimate of the systematic uncertainties, which are frequently 
much larger than the former. In our case, significant systematic er- 
rors arise mostly due to the adoption of an oversimplified mode0. 
The entity of the systematic uncertainties can be quantified, at least 
as a first order approximation, by looking at the discrepancies be- 
tween the recovered parameters for the three data sets best models. 



7 In real datasets of lens galaxies, sources of systematic errors which raise 
particular concern are, for example, the PSF and the sub traction of the 
lens g alaxy surface brightness (see e.g. the detailed study of lMarshall et alj 
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6 DISCUSSION AND CONCLUSIONS 

We have applied cauldron, currently the most advanced code for 
joint gravitational lensing and stellar dynamics analysis of early- 
type galaxies, to a galaxy model with dark matter halo resulting 
from a numerical N-body simulation of galaxy merging. Such a N- 
body system, which we use as lens, significantly violates the two 
major assumptions upon which the algorithm is based, namely axial 
symmetry of the total density distribution and two-integral stellar 
distribution function. The purpose of this crash test is to investi- 
gate how the code will perform or fail in an extreme case, and to 
identify the quantities which can still be reliably recovered. Such 
robust quantities are expected to be recovered with at least compa- 
rable accuracy when cauldron is employed in the analysis of real 
galaxies which, while not necessarily axisymmetric or described 
by a two- (or three-) integral DF, will hardly depart from the code's 
assumptions more severely than the simulated system studied here. 

Further complications can also arise, in general, due to the ef- 
fects of the environment on the lens galaxy. However, the SLACS 
galaxies, despite living in overdense regions as expected for mas- 
sive early-type galaxies, are not found to be significantly affected 
by the contribution of the envir onment or of line -of-sight contami- 
nants (see the in-depth study o f lTreu et al.ll20oi) . Therefore, since 
the environment does not play a major role for lens galaxies at least 
in the redshift range z ~ 0.1 - 0.4, we have not considered this issue 
further in the present paper. 

From the N-body system we have generated three data sets 
corresponding to three orthogonal lines of sight, one of which has 
been chosen to be approximately oriented along the total angular 
momentum of the system, in order to obtain a "face-on" projection 
data set with little or no rotation discernible in the kinematic maps. 
An elliptical Gaussian surface brightness distribution has been con- 
structed and then gravitationally lensed by the mock galaxy in order 
to create the lensing data set. We have also taken into account the 
effect of the PSF and added realistic noise to the simulated data, 
using as a reference the real data set for the SLACS lens galaxy 
SDSSJ232 1-097 studied in C08. 

These data sets have been analysed with cauldron, assum- 
ing as a model an axisymmetric power-law total density distribu- 
tion, with the identical procedure followed in the study of real lens 
galaxies. In the three cases, the recovered best models, obtained 
via maximization of the Bayesian evidence, show clear difficulties 
in reconstructing the observables (in particular the lensed image 
and the velocity dispersion map) up to the noise level. This is a 
consequence of having adopted a very simple model which cannot 
account for the complexity and the lack of symmetry of the true 
density distribution. 

However, the method is still capable of recovering with re- 
markable accuracy several global structural and dynamical charac- 
teristics of the examined system, provided that some information 
about the galaxy rotation is available from the kinematic maps. In 
particular: 

(i) The logarithmic slope y' of the total density distribution is 
a robust quantity which is recovered with remarkable accuracy 
(within less than 10%), even in the case of the face-on projection 
data set. While a power-law model cannot account for a break in 
the actual density profile, indications of its presence will show up 
in the observables as features which are not reproduced by the best 
model (e.g. a much flatter central velocity dispersion). 

(ii) The angle-averaged total density and total mass radial pro- 
files for the edge-on projection best models are found to closely 
follow the corresponding true distributions, within approximately 



an effective radius. Discrepancies are larger for the face-on projec- 
tion case. 

(iii) The best reconstructed inclination angle and flattening of 
the total density distribution have little or no relation with the corre- 
sponding quantities of the N-body system. We conclude that when 
the axisymmetric model is applied to a system with a more com- 
plex geometry, one can not expect to recover reliable information 
about its shape. 

(iv) By adopting the maximum bulge approach (i.e. maximizing 
the contribution of the luminous component, assuming a position- 
independent stellar mass-to-light ratio), it is possible to estimate 
within approximately 10% (in total mass) the dark matter fraction 
of the analysed system, whether the mass ratio is calculated within 
a sphere or within a line-of-sight oriented cylinder of radius = R e . 

(v) When rotation is present in the kinematic maps, global quan- 
tities such as the angular momentum L- and the ratio VI or (a mea- 
sure of ordered vs chaotic motions) recovered from the best model 
describe quite reliably (i.e. within ~ 10% and 25%, respectively) 
the dynamical properties of the system under analysis. The global 
anisotropy parameters /? and y are not correctly estimated, i.e. the 
anisotropy distribution is not robustly recovered by the method un- 
less the analyzed system effectively respects the assumptions of ax- 
isymmetry and two-integral distribution function. The anisotropy 
parameter S, on the contrary, is found to be a robust quantity even 
when such assumptions are violated: we recover the correct value 
within < 15% for the edge-on projection data sets. 

The major conclusion of our study is that the joint lensing and 
dynamics code cauldron can be effectively applied also to the anal- 
ysis of galaxies which deviate from the (quite restrictive) assump- 
tions of axial symmetry and two-integral stellar DF. This result is 
very relevant for the analysis of galaxies at z > 0. 1 for which, due to 
the data limitations, the more powerful methods available at lower 
redshifts are not viable techniques. Several fundamental structural 
and dynamical quantities, in particular the total density slope and 
the dark matter fraction within the region where the data are avail- 
able, can be recovered with good accuracy. Other quantities, such 
as L z and the anisotropy parameter 6, can be reliably recovered pro- 
vided that rotation is detected in the kinematic maps: when this is 
not the case (either because of intrinsic properties of the galaxy or 
because the system is observed face-on) the constraints are much 
looser and one needs to be very sceptical about the outcomes of 
the best model relative to these quantities. We point out that spe- 
cial care should be taken when the best model does not manage 
to reproduce the observables satisfactorily (i.e. at or close to the 
noise level), in particular the lensed image. This is a strong indi- 
cation that the true density distribution deviates significantly from 
the assumptions of the adopted model family, and therefore, while 
there are still reliable recovered parameters (as listed above), more 
delicate quantities such as flattening and inclination angle, as well 
as the reconstructed source, should not be trusted. In these cases, 
the cauldron code can provide a robust but only quite general de- 
scription of the structure of the galaxy, paving the way for the more 
sophisticated modelling techniques which are necessary (together, 
if possible, with a better data set) to achieve a deeper and more 
detailed knowledge of the system. 

Interestingly, large residuals in the lensed images reconstruc- 
tion such as the ones found in this study have never been encoun- 
tered so far in the anal ysis of the SLACS sample of lens galaxies 
jKoopmans et alj|2006l C08) indicating that an axisymmetric sin- 
gle power-law total density distribution constitutes a satisfactory 
model for these systems. 
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